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We present Monte Carlo simulations on Eley-Rideal abstraction reactions of atomic hydrogen 
chemisorbed on graphite. The results are obtained via a hybrid approach using energy barriers de- 
rived from DFT calculations as input to Monte Carlo simulations. By comparing with experimental 
data we discriminate between contributions from different Eley-Rideal mechanisms. A combina- 
tion of two different mechanisms yields, good quantitative and qualitative agreement between the 
experimentally derived and the simulated Eley-Rideal abstraction cross sections and surface config- 
urations. These two mechanisms include a direct Eley-Rideal reaction with fast diffusing H atoms 
and a dimer mediated Eley-Rideal mechanism with increased cross section at low coverage. Such 
a dimer mediated Eley-Rideal mechanism has not previously been proposed and serves as an al- 
ternative explanation to the steering behavior often given as the cause of the coverage dependence 
observed in Eley-Rideal reaction cross sections. 



I. INTRODUCTION 

Molecular hydrogen is the most abundant molecule in 
the interstellar medium (ISM), where it serves as an im- 
portant coolant and as a precursor for the formation of 
more complex molecules. In the cold and dilute inter- 
stellar medium no efhcient gas phase routes exist for the 
formation of II2. Hence, the most likely route is to form 
the molecule on the surfaces of small dust particles. This 
process has been intensely studied over the past years 
m, i, S i, i, i, 0, H, i, m m d- The general con- 
sensus is that molecular hydrogen can be formed by a 
diffusive Langmuir-Hinshelwood mechanism, in the very 
cold regions below ~20 K [13]. At these low tempera- 
tures the hydrogen atoms are physisorbed on the dust 
particles, but can still easily move on the surface. Once 
two atoms meet, they will react to form molecular hydro- 
gen. If the surface temperature becomes too high, the 
residence time of atoms on the surface becomes so short 
that, at the low fluxes in the ISM, the chance of two 
atoms meeting becomes negligible. Observations show 
that molecular hydrogen is also formed in warmer areas 
(> 20 K) like Photon Dominated Regions (PDRs) and 
post-shock regions. A possible mechanism at these con- 
ditions would be through an Eley-Rideal reaction where 
an incoming hydrogen atom reacts with another H atom 
that is chemically bound to the surface. Since a consider- 
able fraction of the interstellar grains is expected to con- 
sist of carbonaceous material, a popular model system for 
this process is hydrogen abstraction from graphite. This 
Eley-Rideal reaction has been studied by Density Func- 
tional Theory (DFT) an d Q uantum Wave packet calcula- 
tions 

[3[il[iS[ii[il, 

lion . Based on these calculations 
several different mechanisms have been proposed to con- 
tribute to the Eley-Rideal abstraction process. These in- 
clude: Direct Eley-Rideal [11, El IH, [S E HI , barrier- 
less abstraction of one hydrogen atom forming part of a 



para-dimcr configuration 20] and abstraction by rapidly 
diffusing H atoms in physisorbed states [2]| . 

However, due to computational limitations such de- 
tailed calculations are unable to take the full complex- 
ity of the H-graphite system into account. Hence, in 
the calculations the flat surface approximation is often 
used and generally either none or only a single C atom 
on the simulated graphite surface were allowed to relax 
during the Eley-Rideal reaction. These simplifications 
are potentially problematic. By not allowing for relax- 
ation of the C atoms on the surface zero barrier reac- 
tion channels, such as sticking of hydrogen atoms into 
specific dimer configuration [2^, [1^, are in some cases 
artificially closed, making it impossible to evaluate the 
contributions from different proposed mechanisms to the 
Eley-Rideal abstraction process. Eg. calculations on the 
Eley-Rideal reaction by Martinazzo and Tantardini fT§] 
show a higher probability for reaction if the incoming H 
atom is 0.5 - 1.5 A away from the target H atom instead 
of a direct hit on. This distance corresponds to the H- 
H distance in an ortho-dimcr configuration (see Fig. 1) 
on the graphite surface (1.42 A). Hence, surface corruga- 
tion and competing sticking reactions must be expected 
to influence the abstraction behavior. In these calcula- 
tions, however, dimer formation was not a possible route 
since the flat surface approximation was used and only 
the carbon atom underneath the original hydrogen atom 
was allowed to relax. Impact parameters corresponding 
to the para-dimer (see Fig. 1) distance (2.84 A) were not 
included in the study. 

Both experimental observations and theoretical calcu- 
lations show that hydrogen atoms preferentially form ad- 
sorbate clusters on the graphite surface. Scanning Tun- 
neling Microscopy (STM) investigations show that at a 
hydrogen atom coverage of ~ 0.5 % more than 75 % 
of all surface configurations are clusters [2^. This indi- 
cates that 85 % of the hydrogen atoms are part of 
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larger clusters. In particular, the existence of hydro- 
gen dimer configurations has been studied theoretically 
[23,[2l0,[2i,01 and experimentally [2i,[23|. Hornekasr 
et al. [24| identified two stable hydrogen dimer configura- 
tions on the graphite surface, an ortho-dimer and a para- 
dimer (see Fig. [T]). DFT calculations show that forma- 
tion of the para-dimer is barrierless and that formation of 
the ortho-dimer has a reduced barrier [25j . which makes 
dimer formation a competing channel to Eley-Rideal ab- 
straction at non-zero impact parameter. 

One set of experiments on the Eley-Rideal abstraction 
reaction forming HD has been reported [1^ [2^ . In these 
experiments a high cross section for the abstraction re- 
action was observed which varied from 17 at low cov- 
erage to 4 at high coverage. This variation in the 
cross section with coverage has been ascribed to a steer- 
ing effect However, other mechanisms could also 
contribute to the high cross section at low coverage. Op- 
tions include the barrier-less abstraction of one of the H 
atoms in a para-dimer configuration [20| and abstraction 
by fast diffusing H atoms [21[ . We propose a further pos- 
sibility, namely abstraction via a hydrogen dimer state, 
where the incoming atom is not immediately thermalized 
and some of its excess initial energy is used to overcome 
the barrier to II2 formation and desorption. The find- 
ing that adsorption into the hydrogen para-dimer state 
is barrierless [25j make this a strong competing reaction 
channel. 

Hence, in total 5 different abstraction mechanisms have 
been proposed to contribute to Eley-Rideal abstraction 
of hydrogen on graphite: 

I. direct Eley-Rideal [Hill. 

II. direct Eley-Rideal with steering [TgI. [TtI. [2§| . 

III. preferred direct Eley-Rideal with hydrogen atoms 
part of a para-dimer . 

IV. a dimer mediated reaction where the incoming 
atom is first adsorbed into a dimer configuration 
and, before thermalizing to the substrate temper- 
ature, reacts with the other atom in the dimer to 
form H2 and desorb. 

V. Eley-Rideal reactions by fast diffusing H atoms in 
the physisorption state [2l|. 

However, due to the simplifications needed in the com- 
plex DFT and quantum wave packet calculations a quan- 
titative comparison between experimental data and the- 
ory has not been possible. In this paper we employ a 
hybrid approach in which we include the findings of ab- 
initio DFT calculations in a Monte Carlo simulation pro- 
gram and then simulate Eley-Rideal abstraction experi- 
ments on a more realistic surface area with more complex 
hydrogen adsorbate configurations than what is possible 
in the DFT and wave packet calculations. Through this 
quantitative approach we aim to discriminate between 
the contributions of the different underlying mechanisms 
for hydrogen abstraction. 
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FIG. 1: The two stable dimer configurations found by 
Hornekffir et al. 12411. 
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FIG. 2; The adsorption sites used in the Monte Carlo simu- 
lations. 



II. MONTE CARLO MODEL 

The Monte Carlo simulation program is a so-called 
lattice-gas simulation program, where the atoms are con- 
fined to an adaptive grid depending on the initial cover- 
age. The hydrogen atoms can chemisorb to the graphite 
surface at the sites directly on top of the carbon atoms 
as indicated in Figure [21 The puckering of the carbon 
atom upon chemisorption is included indirectly by using 
the barriers for chemisorption and binding energies from 
DFT calculations that allow this motion. Interaction be- 
tween two hydrogen atoms is accounted for in a similar 
way as will be discussed in the following sections. The hy- 
drogen atoms can physisorb to the surface without a bar- 
rier. Since the potential energy surface for physisorption 
is rather fiat, there are probably no specific physisorp- 
tion site, however since the lattice-gas model forces us to 
choose specific sites, we confine physisorption to the sites 
directly above the carbon atoms and an additional site 
at the center of the ring (Figure [2]). 

The simulation starts with a clean graphite surface. 
The first event will be a deposition attempt of the first 
atom. The time at which this attempt will occur is 



^dcp 



In (X) a 



+ t 



(1) 



where a is the density of sites, X is a random number 
between and 1, / is the hydrogen fiux in atoms per 
time per area, and t is the current time. Deposition 
times for subsequent sticking events are determined by 
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the same expression (Eq. [T]) . How the atoms will bind 
to the surface is determined by another random number 
and depends on the barrier for chemisorption in that spe- 
cific position. If the site is already occupied by another 
hydrogen atom, an abstraction reaction is considered by 
comparing a Boltzmann factor including the abstraction 
barrier against a random number. Upon reaction both 
atoms will leave the surface in the form of H2, else the 
incoming atom is deflected. The abstraction channel can 
be closed by making the barrier infinitely large. 

Once hydrogen atoms populate the surface, they can 
diffuse, desorb, or recombine with other atoms to form 
H2. For each surface hydrogen the time at which they 
will undergo one of these events is determined by 



ln(X) 



R 



dif 



+ t 



(2) 



with R 



■dif 



the rate for diffusion, desorp- 
tion, and reaction of atom i, respectively. Another ran- 
dom number determines which of the three events occur 
according to their relative probability of occurrence. Re- 
action between two chemisorbed atoms can only occur if 
they form a dimer configuration. The rates are given by 



R = V cxp — — — 



E 
kT' 



(3) 



where v is the attempt frequency which is assumed 
to be 10^^ Hz for physisorbed atoms and 10^^ Hz for 
chemisorbed atoms [2J| and T' is the 'temperature' of 
the atom that is involved. 

In the experiments that we aim to reproduce, the 
atoms arrive at the surface at normal incidence and 
at a temperature around 2000 K which is much higher 
than the surface temperature. Furthermore, the atoms 
will gain energy due to the high binding energy if they 
chemisorb. The atoms will not be thermalized instanta- 
neously, but will most likely gradually lose their energy 
to the substrate. The dissipation of the excess energy 
into the substrate is expected to be exponential [s^l- To 
reduce the computation complexity of the model we em- 
ulate this exponential energy loss by a simpler expression 
and use the following function to describe the tempera- 
ture 



T'(t) 



{1 + B{t-QY 



(4) 



where Ts and Tstart are the surface and starting temper- 
ature respectively, is the time at which the atom has 
adsorbed on the surface. The sensitivity of the model to 
the exact functional form was checked and found to be 
small. For most cases we use Tgtart = 2000 K, but also 
higher values of 7000 K and 10,000 K are considered, 
matching the binding energy of a monomer. The param- 
eter B can be chosen freely. Fig. [3] studies the effect of 
this parameter. It displays the result of series of 10,000 
simulations of one deposition event for a particular value 
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FIG. 3: Influence of the B parameter on the sticking fraction 
(top) and the relaxation time (bottom) of the hot atom. For 
detailed explanation see text. 



of B. The top panel indicates the percentage of depo- 
sition attempts that resulted in sticking. The sticking 
fraction slowly approaches the monomer sticking barrier 
probability of 41 % for increasing B. A value of B above 
^ 10^ s^^ for Tstart = 2000 K is needed to get an initial 
sticking co-efficient in agreement with the experimental 
findings |28i|. Finally, the bottom panel gives the relax- 
ation time. For comparison Shalashilin and Jackson (soj 
found for a hydrogen atom on a Cu(lll) surface that the 
thermal relaxation time is around 4 ps. This corresponds 
to a B of 8 X 10^° s-i for Tstart = 2000 K and 2 x 10" 
s-i for Tstart = 10, 000 K. 

In the brief moment after a deposition, the rates in 
Eq. [3] are time dependent due to the decreasing tem- 
perature and Eq. [5] cannot be used. Instead we use 
the method by Jansen et al. [3l| to determine i*. This 
method makes use of 

-ln(X)= RUt)dt+ f i?J,es(t)dt-f /* i?J.ecWdt, 
Jt Jt Jt 

(5) 

that transforms to Eq. [2]if all rates are time dependent. 
To obtain ti the expression has to be solved. Using 



i/exp — ]dt (6) 

t \ start 



" ^"^^*-'-erf (l + B(t-U)J-^\7) 



2B\ E 



Ts, 



with erf the error function, Eq. [5] becomes 

- In (X) = 17dif (f ) + l^dcs {t' ) + rircc it' ) . (8) 
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This can be solved numerically to ti using the Newton- 
Raphson method IS^j], since both O and ^ decrease 
monotonically. Notice that different functions for T' will 
result in a different expression for ri(i*). As the order in 
t increases, solving r2(i') becomes more computationally 
expensive. 



TABLE I: The sticking barrier for different configurations 



Configuration 


Bstick [mcV] 


para-dimer 


0.0 


ortho-dimer 


0.1 


other chemisorption 


site 0.15 



physisorption 0.0 



III. INCLUDED REACTIONS AND ENERGY 
BARRIERS 

The previous section described the general Monte 
Carlo algorithm. For all processes energy barriers are 
needed to determine the corresponding transition proba- 
bilities. These barriers are taken from independent DFT 
calculations ^33\ of the binding energies of the different 
configurations that are considered, sticking trajectories 
and diffusion trajectories. Since many possible config- 
urations are formed during the simulations, especially 
for high coverages, including all these different possibili- 
ties explicitly would make the Monte Carlo program very 
slow and would require a huge set of barriers that all have 
to be calculated independently. To overcome this prob- 
lem a number of simplifications are introduced resulting 
in the following sets of energy barriers: 



A. Physisorbed atoms 

For the physisorption binding energy, the value of ~ 40 
meV is used based on results from selective adsorption 
experiments [s^l- For diffusion an activation energy of 
4 meV taken from [2l| was used. We will come back to 
this diffusion rate at the discussion of mechanism IV. For 
physisorbed atoms T'(t) = Tg was used. 



B. Sticking 

Numerous different barriers for sticking of an isolated 
hydrogen atom (a monomer) into the chemisorption site 
on the graphite surface have been given. Jeloaica and 
Sidis found a barrier of ~0.2 eV using the coronene 
molecule as a model of a graphite surface Sha 
et al. obtained a barrier slightly above 0.2 eV using a 
slab super cell with 4 layers each containing 8 carbon 
atoms [IB, [13] ■ Hornekaer et al. [2^ found a barrier for 
chemisorption into a monomer of 0.15 eV using a single 
layer super cell containing 32 carbon atoms. The influ- 
ence of adding a second carbon layer was investigated and 
found to be negligible. Using the same model surface a 
sticking barrier of 0.1 eV into the ortho-dimer and eV 
into the para-dimer configurations were found. All possi- 
bilities for sticking into a trimer state, starting from the 
para-dimer, were seen to have non-zero sticking barriers 
between 0.1 and 0.15 eV. Adding a fourth atom resulting 
in a triple para-dimer configuration again did not exhibit 



TABLE II: The total binding energy for different configura- 
tions 



Configuration^ 


Ebind [meV] 


monomer 


-0.8 


para-dimer (G) 


-2.9 


ortho-dimer (A) 


-2.8 


meta-dimer (E) 


-0.8 


trimer (A, G) 


-3.9 


trimer (G, I) 


-3.7 


trimer (D, G) 


-3.5 


trimer (B, G) 


-3.5 


trimer (F, G) 


-3.3 


tetramer (A, G, H) 


-5.9 


tetramer (A, B, E) 


-5.9 


tetramer (A, D, G) 


-5.7 


tetramer (A, B, D) 


-5.5 


tetramer (B, E, G) 


-5.1 


tetramer (A, B, C) 


-5.0 


tetramer (A, B, G) 


-4.8 


tetramer (A, B, F) 


-4.5 


1 See Fig. 


m 



a barrier [25|. Barrierless sticking into the para-dimer 
state was also found by Rougeau et al. (26j . 

Even though there is a clear dependence on the local 
configuration of the impact site, we decided only to in- 
clude variations in sticking barriers for dimers in the sim- 
ulation. The program determines if the incoming atoms 
can form a dimer ignoring the larger configuration it 
might be part of and determines the barrier for sticking 
accordingly. This assumption will cause deviations be- 
tween simulations and experimental results at high cov- 
erage, since it overestimates the formation of trimer con- 
figurations that contain para-dimers. No sticking barrier 
was used for physisorption of H atoms. The barriers for 
sticking used in the simulation are summarized in Table 
III 

Section tests the different mechanisms to their sen- 
sitivity to various input parameters. The monomer stick- 
ing barrier is one of them. 



C. Binding energies 

The binding energies for different configurations of 
chemisorbed H atoms are displayed in table II. The dif- 
ferent configurations are displayed in Fig. [J] These values 
are based on DFT calculations reported in [2^ [s^ . 
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the same relation. A least-squares fit resulted in 

Earn = 0.5AE + 1.04 eV. (9) 

In order to obey detailed balance, or microscopic re- 
versibility, transition probabilities for diffusion should 
fulfill 



P. 



exp 



kT 



(10) 



Using Eq. [31 it can be shown that the empirically found 
relation, Eq. [HI follows this requirement. 



E. Desorption 



FIG. 4: Schematic guide to obtain the configurations used in 
Table [TT] A configuration is made up from hydrogen atoms 
positioned on top of the black carbon atom and the atoms 
indicated by the characters given in Table [m In this way, an 
ortho-dimer is represented by the character (A). 




O dimers 
-I- tetramers 
— 0.5 AE + 1.04 



AE (eV) 



1 



FIG. 5: The diffusion barrier as a function of the energy differ- 
ence between the initial and final configuration. The diffusion 
barriers between dimer and the tetramer configurations follow 
the same linear dependences. 



Finding a general expression for the desorption en- 
ergy of the individual atoms from a configuration is less 
straightforward. The binding energies are determined for 
the configuration and not for the individual atoms. Again 
calculating all possible desorption pathways would not 
be feasible. We base the desorption on the total bind- 
ing energy of the configuration of n atoms, ii^bind- The 
desorption energy is then 



P _ -^bind j-1 
-C'dcs — r -CVstick 



n 



(11) 



where -Estick is the barrier for sticking in the same po- 
sition, -©stick and -Ebind can be found in Tables H] and 
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F. Thermally activated H2 formation 

Formation of H2 via reaction between two chemisorbed 
hydrogen atoms to gaseous molecular hydrogen occurs 
from the dimer states and has to overcome barriers of 
2.49 eV for the ortho-dimer and 1.4 eV for the para- 
dimer state [ISl- Reaction barriers from the isolated 
dimer states are used for all configurations. This approx- 
imation will again lead to inaccuracies for simulations at 
high coverage. 



D. Diffusion 

For several dimer and tetramer configurations the in- 
dividual binding energies and the transition barriers be- 
tween some of the configurations were determined via 
DFT calculations. The binding energy of the individual 
dimer and tetramer configurations are given in Table [Til 
The tetramer binding values and diffusion barriers are 
taken from [s^ . Figure [5l plots the diffusion barrier as 
a function of the energy difference of the configurations. 
As the figure clearly shows there is a linear relation be- 
tween the two quantities and dimers and tetramers follow 



G. Abstraction 

Five different Eley-Rideal abstraction mechanisms are 
considered in the model: 

I. direct Eley-Rideal with all hydrogen atoms regard- 
less of their local configuration. 

II. a simple version of Eley-Rideal with steering where 
a direct Eley-Rideal reaction is allowed not just for 
H atoms impinging on an already occupied site but 
also for H atoms impinging on adjacent sites. 

III. direct Eley-Rideal with only hydrogen atoms part 
of a para-dimer. 
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FIG. 6: The amount of HD molecules per time leaving the surface in ML/s as a function of time for different pre-exposures 
of D. At t = the hydrogen beam is switched on. For (I) mechanism L Direct Eley-Rideal (BER.mono ~ i5ER,dimer = 9 meV), 
(II) mechanism II: Direct Eley-Rideal with steering (i?ER,mono ~ -EER.dimcr ~ 9 meV, s — 1), (III) mechanism III: Dimer 
Eley-Rideal {Eer, .mono — CXD; -EERjdimcr = meV), (IV) mechanism IV: Dimer mediated Eley-Rideal (B = 1 x 10^ s ^), and (V) 
mechanism V: Direct Eley-Rideal with fast diffusion of physisorbed atoms (E'er mono = £er dimer ~ 9 meV, Rdic — 5 x 10^^ 
s-i). 



IV. a dimer mediated reaction where the incoming 
atom is first adsorbed into a dimer configuration 
and, before thermalizing to the substrate tempera- 
ture, forms H2 and desorbs. 

V. direct Eley-Rideal together with a high diffusion 
rate of the atoms in the physisorption state 

The influence of the different Eley-Rideal mechanisms 
can be controlled by using different barriers for the direct 
Eley-Rideal reaction and different values of the thermal- 
ization parameter B. Furthermore, a simple version of 
steering is implemented with an adjustable parameter s 
as will be described below. The values of the parameters 



used for simulating the different mechanisms are given in 
Table Uni 

The height of the abstraction barrier is somewhat un- 
certain due to the limitations in the accuracy of DFT 
calculations. Morisset et al. [18] found a barrier for di- 
rect abstraction of a monomer just below 10 meV. Others 
also found low barriers to abstraction J/^, 19]. Follow- 
ing Morisset et al. we employ a barrier of 9 meV but 
later investigate the effect of changing the value of the 
barrier. The same group also found that the abstraction 
reaction with one of the hydrogen atoms in a para-dimer 
can proceed without barrier |20| . 
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TABLE III: Model parameters for the different abstraction mechanisms 



Mechanism 


-EER.mono [meV] 


-EEH,,dimcr [mcV] 


B [s-^] 


Tstart [K] 


s 


Rdii[s''] 


I (Fig. HI) 


9 


9 




2000 





7.2 X 10'^ 


II (Fig. nil) 


9 


9 




2000 


1 


7.2 X 10^^ 


III (Fig. niii) 


oo 







2000 





7.2 X 10^2 


IV (Fig. HIV) 


(X 


oo 


10^ 


2000 





7.2 X 10^^ 


V (Fig. HV) 


9 


9 


f0^2 


2000 





5.0 X 10^^ 



IV. COMPARISON WITH EXPERIMENTS 

The Monte Carlo simulation results are compared with 
the experimental data on the Eley-Rideal reaction pre- 
sented by Zecho et al. [11] . These experiments consist of 
two phases. First, a graphite substrate was exposed to a 
normal incidence deuterium beam with a flux of 3.8 x 10^^ 
atoms cm~^s~^ at 150 K. The exposure time was varied 
to give a range of initial coverages. The temperature of 
the atom beam was 2000 K. During the second phase the 
pre-exposed substrate was exposed to a normal incidence 
hydrogen beam and the formed HD molecules were mea- 
sured using a mass spectrometer. The cross section of the 
deuterium abstraction with hydrogen, a, was determined 
from these spectra using 

= a$[Dad,o] exp (-a$t) , (12) 

with $ the H flux and [Dad.o] the initial D coverage. For 
the derivation of this equation, we refer to [1^. The 
expression is obtained assuming only direct Eley-Rideal. 
Since the reaction mechanism is implicitly included in 
the cross section, it is hard to directly interpret the re- 
sults obtained in this way, but we will use the method 
to compare the experimental data with our simulation 
results. 

The hydrogen adsorbate configurations found in the 
Monte Carlo simulations are also compared to the ad- 
sorbate configurations observed in STM experiments. In 
particular the fraction of hydrogen atoms in dimer config- 
urations or larger clusters in simulation and experiment 
are compared. 



V. RESULTS 

Several simulation runs were performed using a sim- 
ilar set of pre-exposures as used by Zecho et al. [1^. 
The H and D fluxes were 4.8 x 10^^ cm~^s~^ at a tem- 
perature of 150 K. 250 x 125, 500 x 250, and 750 x 500 
chemisorption sites were used depending on the initial 
D coverage. This corresponds to an array of 250 x 250, 
500 X 500, and 750 x 750 to accommodate the extra ph- 
ysisorption sites. If the noise was still considerable for 
arrays of 750 x 500 multiple runs with different random 
seeds were made. In order to test to what extent the 
five proposed abstraction mechanisms contribute to the 
Eley-Rideal abstraction process comparisons to the ex- 



perimental results, with different parameter sets (see Ta- 
ble |Tn| were performed. 

Figure [6] shows d[HD]g/dt for the five different abstrac- 
tion mechanisms. Panel (I) only includes the direct Eley- 
Rideal with a barrier of 9 meV independent of the con- 
figuration on the surface, panel (II) investigates the in- 
fiuence of steering, panel (III) only allows an Eley-Rideal 
reaction, with a zero barrier, if the surface atom is part 
of a para-dimer configuration, panel (IV) uses the slower 
thermalization with _B = 1 x 10^ s~^ to test the dimer 
mediated mechanism, and panel (V) uses very fast dif- 
fusion. As d[IID]g/d< is plotted on a logarithmic scale, 
the curves should be linear according to Eq. [T^ with the 
slope the cross section times the fiux. The flux has the 
same value throughout the simulations. 



The key results of the simulations are summarized in 
Table IIVI and compared to the experimental findings. 
Each entry consists of nine individual simulations. The 
top twelve rows represent the runs that allow only a single 
mechanism. The second section summarizes simulation 
series of combinations of mechanisms and the third part 
gives the results with a different sticking barrier and will 
be discussed in Section IVll The final entries give the ex- 
perimental results for comparison. The difference in cross 
section between high and low coverage, the saturation 
coverage, and the linearity of the initial signal in initial 
coverage are indicators of the agreement with the exper- 
iment. The latter linearity is measured by the Pearson 
correlation coefficient that gives one for perfect correla- 
tion and zero for no correlation. For a better comparison 
between the different simulation runs all cross sections at 
low coverage are given for ~ 0.013 ML which is generally 
achieved after an exposure of 0.03 ML (the lowest con- 
sidered pre-exposure). Simulations including mechanism 
IV (dimer mediated reaction) need a longer exposure to 
reach 0.013 ML whereas simulations with mechanism V 
(fast diffusion) reach it faster than 0.03 ML. The cross 
sections at high coverage are taken at the highest consid- 
ered pre-exposure of 2.9 ML. The tenth column gives the 
corresponding coverages. In almost all cases the satura- 
tion coverage has been reached at this point. The last 
column of the table indicates the fraction of atoms that 
is part of a dimer configuration or larger cluster at a cov- 
erage of 0.5 % at conditions similar to the ones used by 
Hornekaer at al. [23 |. 
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Mechanism I: Direct Eley-Rideal 

The five plots in Fig. [S]all show a very different behav- 
ior. Panel (1), which tests the direct Eley-Rideal mecha- 
nism, shows a linear relation between In (d[HD]g/dt) and 
the deposition time over the whole time range. The slope 
appears to have some dependence of the initial D cover- 
age. This would result in a cross section that is dependent 
of initial coverage. Figure [7] plots, among other quanti- 
ties, the cross section which indeed shows some coverage 
dependence. It is however much less than seen in the 
experiments [2^. The cross section and the value of the 
initial signal is obtained by fitting the curves shown in 
Figure to Eq. [T^ Since for some parameter choices 
only the first section is linear, only these points are in- 
cluded in the fit. The number of points included in the fit 
can have a strong effect on the final result, especially if 
the linear part is very short. This is reflected by the error 
bars. Various runs with different initial seed showed that 
this error should at least be 0.5 A^. The initial signal de- 
pends linearly on the deuterium coverage, which agrees 
with Figure 3c in [2^ . Figure [7] also displays the total 
HD yield. Since this value is comparable with the initial 
D coverage it is clear that only a negligible fraction of D 
atoms remains on the surface after hydrogen exposure. 
This is in good agreement with the experimental find- 
ings f28| . The dimer fraction is very low in contrast with 
the experiments as can be seen in the first row in Table 

CYl 



Mechanism II: Eley-Rideal with steering 

A simple version of Eley-Rideal abstraction with steer- 
ing is implemented by allowing atoms that land on empty 
sites close to chemisorbed hydrogen atoms to attempt a 
reaction with these atoms. For the reaction barrier the 
same i^ER.mono is used as for the direct hit. Eley-Rideal 
via a direct hit (mechanism I) is also allowed. If an atom 
lands on an empty site, the possibility of steering is con- 
sidered by comparing a random number between and 
1 against parameter s. If the random number is smaller 
than s, then steering is activated. In the case of steer- 
ing, the three neighbouring sites are checked for react- 
ing species. If two or more different atoms are found, 
another random number determines with which of the 
atoms the incoming H atom will undergo the reaction at- 
tempt. Simulations for s = 1 and 0.5 were performed. 
The results of these two simulation runs are summarized 
in Fig. [5] and the rows indicated by II in Table IIVI and 
the s = 1 simulation run is shown in Figure [6]-II. Fig. [8] 
shows the initial coverage and yield after 200 seconds of 
H exposure as function of the initial D pre-exposure. No- 
tice that the yield is given in ML in these graphs and not 
as a percentage. It further gives the cross section and 
initial signal as a function of the initial coverage. For 
s = 1, a reasonably high cross section at low coverage 
is obtained (14.7 ± 1.4 A^), but the saturation coverage 
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is very low (0.10 ML) and also the cross section at this 
coverage is too high. The series which simulates steering 
with s — 0.5 gives somewhat lower cross section result- 
ing in better agreement with the experimental findings at 
high coverage and too low a cross section at low coverage. 
Furthermore, an almost linear decrease in cross section 
is found as a function of coverage in disagreement with 
the experimental findings. Again the dimer ratio is too 
low. 



Mechanism III: Dimer Eley-Rideal 

The curves in panel (III) are not linear and show a 
maximum at later times for the low coverages. Since 
this panel presents the results with the dimer Eley-Rideal 
mechanism, this maximum indicates the time at which a 
maximum of D containing dimers is formed. This be- 
havior contradicts the experiments [1^ where the [HD] 
signal decreases over time and is therefore not likely to 
be the primary mechanism involved in the abstraction. 
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Mechanism IV: Dimer mediated abstraction 

Panel (IV) presents the results for the dimer mediated 
mechanism. Here a clear linear dependence can be ob- 
served for very early times with a decreasing slope for 
increasing initial coverage. To obtain this graph large 
arrays upto 1500 x 3000 sites were used due to the low 
sticking rate at low coverage. Notice that the y-axis range 
of this panel is different from the others. Figure [5] studies 
this set of simulations more closely. The two upper panels 
show that not all D is converted into HD, but that a sub- 
stantial amount of deuterium still resides at the surface 
after the 200 seconds exposure. This is in contrast with 
the experimental finding that all deuterium is abstracted 
after a 3 ML dose. Figure [H] also shows the cross section 
determined from Fig. [6]-II. Again the bars indicate the 
uncertainties in obtaining the slope. Since the linearity 
is much less than for mechanism I, the error bars are 
larger than in Fig. [T] In contrast with this latter figure, 
Figure [9] shows a strong dependence of the cross section 
on the initial deuterium coverage. This corresponds with 
the trends observed by Zecho et al. HI] . The values are 
also in the correct range, although the cross section is too 



high at low coverage and too low for high coverage. Both 
the cross section and the initial signal indicate that the 
dimer mediated mechanism becomes very inefficient for 
high coverages. The initial signal even decreases for in- 
creasing surface coverage indicating that at high coverage 
it is hard to make dimers due to the decreasing number 
of available sites. As the coverage increases, either by H 
or by D atoms, the mechanism becomes less efficient and 
the deuterium atoms will not be abstracted. This also 
results in a very high saturation coverage. 

The dimer fraction is higher as compared to the other 
mechanisms, although still much too low. The increase 
in dimer fraction is due to the fact that many of the 
chemisorbed atoms desorb due to the initial energy. Since 
atoms in dimer position are less likely to desorb because 
of their higher binding energy, this results in an elevated 
dimer ratio. 



Mechanism V: Fast diffusion 

As discussed above, the barriers to diffusion in the 
chemisorbed state are quite high making this process es- 
sentially negligible at least at low coverage. Diffusion in 
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the physisorbed state is, however, a completely different 
matter. If a simple expression for an activated process is 
assumed, then an atom thermalized to a surface temper- 
ature of 150 K will make around 20 hops, which corre- 
sponds to a distance of 5 A before desorbing. However, 
as Bonfanti et al. [2l| pointed out diffusion is not an ac- 
tivated process, since the activation barrier of 4 meV is 
negligible at these temperatures. This asks for a differ- 
ent estimation of the moving rate. As a first approach 
we took the diffusion coefficient obtained by Bonfanti ct 
al. [12]. This results in a moving rate of 1.3 x 10^"^ s^^. 
The upper limit for diffusion will be free movement in two 
dimensions according to the average gas phase velocity. 
This results in a diffusion rate of 5.0 x 10^'^ s~^ at 150 K. 
We will perform simulations using both rates. The results 
for the highest rate and i?ER,mono = 9 meV are shown in 
Figures [6]-V and [TOl Results of other simulation series 
with different combinations of i?ER,.mono and the diffu- 
sion rate are summarized in Table [TVl A clear increase 
in the dimer fraction can be observed for increasing dif- 
fusion rate, although still slightly too low compared with 
experiments. The cross section curve as a function of the 
coverage has the correct dependence although it is still 



not high enough at low coverage. Furthermore the initial 
signal as a function of coverage does not have a linear 
dependence but exhibits a kink around 10 %. At this 
coverage the number of trimers starts declining in favour 
of the number of tetramers. Most of the configurations 
(60 %) are now a configuration of three or higher. This 
non-linearity could well be an artifact of the simplified 
treatment of reaction and sticking for these configura- 
tions which makes the model inaccurate at high coverage 
as discussed above. 



Combined mechanisms 

The experiments cannot be explained by one of the 
five mechanisms separately. Mechanism I is not efficient 
enough at low coverages, mechanism II leads to very low 
saturation coverages, mechanism III shows a maximum 
at t > 0, mechanism IV is not efficient at high coverages. 
Moreover, all these four mechanisms show very low dimer 
fractions compared with experiments. Mechanism V has 
an increased number of dimers and a higher cross section, 
but again this is not high enough to match the experi- 
mental observations. We therefore investigate different 
combinations of the five mechanisms. 

Again the results are summarized in Table ITVl The ta- 
ble clearly shows that a combination of mechanisms IV 
and V gives the best results, both in terms of dimer ra- 
tios and cross section. The best agreement is obtained 
for the high diffusion rate and i?ER.mono = 9 meV (in 
bold face) . The detailed analysis of these simulations is 
given in Figure [TlJ The cross section has the correct cov- 
erage dependence, both in trend and in value. Also the 
saturation coverage of 0.4 ±0.2 ML found experimentally 
is reproduced. The major discrepancy is again the initial 
HD signal which does not dependent linearly on the ini- 
tial cross section, but as mentioned earlier this can be an 
artifact due to the simplified assumptions made for stick- 
ing and reaction for complex configurations, which make 
high coverage simulations inaccurate. The dimer ratio is 
still below the experimental values, but much higher than 
all previous series. The dimer ratio includes both ortho 
and para dimers, but the main contribution is from the 
para dimers (40 times more at 0.1 %), since these sites 
have no barrier for sticking. The residence time of the 
physisorbed atoms at ortho sites is not long enough for 
the atoms to easily enter the chemisorbed state. As the 
dimers have been observed in more equal ratios, the rate 
for crossing this barrier is probably higher, possibly due 
to a lower barrier or higher pre-exponential factor than 
used in the present model or due to a longer residence 
time at the ortho site because of the irregularity in the 
lattice caused by the first hydrogen atom. The reason 
that this combination of mechanisms give higher cross 
section at low coverage is because more reaction events 
go via a dimer state. In this state an atoms gains some 
energy due to the strong binding energy and this extra 
binding energy is used to overcome the barrier for reac- 
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tion. 

Other combinations show less agreement with the ex- 
periments on at least two points of comparison. 



VI. SENSITIVITY OF THE MODEL ON 
PARAMETERS 

By comparing the results listed in Table IIVI it is clear 
that a combination of mechanisms IV and V reproduces 
the experimental findings quite well. This section will 
discuss the influence of the parameter settings on the 
final result. The table already showed four different pa- 
rameter choices varying the diffusion rate and the barrier 
for abstraction. The diffusion rate has, as mentioned be- 
fore, a clear influence on the number of formed dimers, 
but it appears to have very little effect on the obtained 
cross sections for the mechanisms IV and V. Lowering 
it to a very low value of 7.2 • 10^^ s^^, which turns it 
into a combination of mechanisms I and IV, results in 
a lower cross section at low coverage. We found that 
reducing the diffusion rate even further eventually leads 
to a completely flat cross section and a very high Pear- 
son correlation coefficient. Also for a pure mechanism V 
a higher diffusion rate gives a stronger coverage depen- 
dence of the cross section as was suggested by Bonfanti 
et al. 

The abstraction barrier has no effect on the dimer ra- 
tio, but it influences the cross section at low coverage. It 
generally gives a lower cross section for a higher barrier. 
For a pure mechanism I also a barrierless abstraction was 
considered (see Table ITVT) . The cross section at low cov- 
erage did not raise to the values found in the experiments 
in this case. 

The parameter B, which controls the thermalization of 
the sticking atoms, has a very narrow parameter range 
if only mechanism IV is considered. Lower values than 
1 • 10^ s~^ result in very low sticking at low coverage 
and higher values will have only a limited effect on the 
cross section. If mechanism IV is used in combination 
with mechanism V the sticking will be increased because 
of the diffusion of the physisorbed atoms and also lower 
values than 1 • 10® s~^ can be used leading to a higher 
cross section for low deuterium exposure and high dimer 
ratios. The linearity in the initial HD mass spectrometer 
signal is however decreased. 

Adding mechanism III as a third mechanism in the sim- 
ulations had a negligible effect on the results. The cross 
section and initial signal curves remained unchanged 
within their uncertainties. 

The parameter, Tgtart, which controls the initial energy 
of hydrogen atoms, has a different effect on the final re- 
sults for the pure mechanism IV and the combination of 
IV and V. This parameter cannot be controlled indepen- 
dently of B, since the sticking probability is constrained 
within a certain window (see Fig. [3]). For a pure mech- 
anism IV a increased value of Tgtart results in a cross 
section which is too low and only weakly dependent on 



coverage. The correlation coefficient for the initial HD 
signal is higher as compared to Tstart = 2000 K. This 
parameter combination further shows an extremely high 
dimer ratio. Only combinations of mechanism were found 
to yield similar high values. It appear that the high ini- 
tial energy is only used to remove monomers and not 
for dimers to react or to desorb. The relaxation time 
is probably to fast for these events to occur. If a more 
elaborate implementation of this mechanism was consid- 
ered, where the initial energy of atoms in dimers is high, 
because of their higher binding energy, the reaction from 
dimers would go up and the dimer ratio go down. The 
results would then probably closer resemble the Tgtart = 
2000 K results. For a combination of mechanisms Tgtart 
has less of an effect. Here only the saturation coverage 
shows a large change. 

Finally we checked the influence of increasing the 
monomer sticking barrier to the value of 0.25 eV which is 
an upper value for the sticking barrier reported by other 
authors [l^ [13, [11] ■ This results in a higher cross sec- 
tion at low coverage and a higher dimer ratio in all sim- 
ulated cases. However, the correlation coefficient of the 
initial mass spectrometer signal as a function of coverage 
goes down. In general this increase in the sticking bar- 
rier results in a better agreement with the experiments. 
Although the agreement of the pure mechanism V is im- 
proved much, the combination of mechanism IV and V 
remains the best option for both low and high values of 



VII. CONCLUSION 

The Eley-Rideal abstraction of atomic hydrogen 
chemisorbed on the graphite surface has been studied 
via a hybrid approach using energy barriers derived from 
DFT calculations as input to Monte Carlo simulations. 
Through comparison with experimental data we discrim- 
inate between the contributions from different proposed 
Eley-Rideal mechanisms. Good quantitative and qual- 
itative agreement between the experimentally derived 
and the simulated Eley-Rideal abstraction cross sections 
are found if two different Eley-Rideal abstraction mecha- 
nisms are included. One is a direct Eley-Rideal reaction 
with very fast diffusion of physisorbed H atoms leading to 
the formation of hydrogen dimer configurations, while the 
other is a dimer mediated Eley-Rideal mechanism with 
increased cross section at low coverage. Such a dimer me- 
diated Eley-Rideal mechanism has not previously been 
proposed. 

The effect on abstraction of fast diffusing physisorbed 
H atoms was first considered by Bonfanti et al. ^2ll| , who 
suggested that diffusion of the physisorbed atoms could 
explain the high coverage dependence of the cross sec- 
tion. We tested this mechanism and it was indeed found 
that diffusion plays an important role, not only for the 
abstraction reaction but also to explain the high occur- 
rence of dimers found by STM measurements. However, 
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TABLE IV: Simulation results for different parameter settings 
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very little is known about the exact energetic landscape 
in the vicinity of a chemisorbed hydrogen atom, espe- 
cially when the second atom is only weakly physisorbed. 
Since these weak interactions appear to be very impor- 
tant, further study would be desirable. 

Furthermore, we investigated the effect of steering 
[H, [l3| as a possibly alternative route to reproduce the 
experimental findings. The results presented here show 
that a simple steering mechanism which just results in 
increased cross sections for abstraction reactions is not 



sufficient to reproduce the experimental results. This 
finding does, however, not rule out steering as an impor- 
tant mechanism. It has for example been suggested that 
a coverage dependent steering effect could result from 
the creation of and interaction with electron-hole pairs 
on the surface [s^. Further investigations of steering 
should, however, consider the presence of preferred bind- 
ing sites in the vicinity of adsorbed hydrogen atoms. In 
the present case, sticking in a dimer position with low or 
no barrier is a competing reaction to the steered Eley- 
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Rideal reaction. Hence, the described dimer mediated 
Eley-Rideal mechanism offers a competing mechanism to 
steering. Full quantum mechanical calculations where 
both channels, dimer formation and hydrogen abstrac- 
tion, are accessible, are needed to settle this case. 
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